xxxxxxxxxx9
1
begin2
using BenchmarkTools3
using CSV4
using LaTeXStrings5
using NearestNeighbors6
using Plots7
using PlutoUI8
using Random9
endxxxxxxxxxx1
1
br = html"<br>";xxxxxxxxxx1
1
bigbr = html"<br><br><br><br><br>";xxxxxxxxxx1
1
html"<button onclick=present()>Present</button>"Modeling and Optimization
An epidemic example
Want all the details?
This talk is based on the free, online course "Computational Thinking with Julia" available from juliaacademy.com
Lecture notes and problem sets (.ipynb and .pdf) are available at https://github.com/mitmath/6S083/
An updated and expanded version of the course called "Introduction to Computational Thinking" is available at https://computationalthinking.mit.edu/Spring21/
Further reading:
Algorithms for Optimization by Kochenderfer & Wheeler
Algorithms for Decision Making by Kochenderfer, Wheeler, & Wray
The Procedure
1. Develop a computer model to simulate reality
2. Establish a theoretical foundation to explain results
3. Compare the computer model to the theory
The Random Walk Model
Let N "agents" randomly move within a square of side L.
For true "agent-based modeling" (ABM), use Agents.jl.
First define the walkers and their methods.
Use an "enumerated type" to keep track of infection status
Susceptible -> Infectious -> Recovered
xxxxxxxxxx1
1
InfectionStatus S I RR::InfectionStatus = 2xxxxxxxxxx1
1
Rxxxxxxxxxx15
1
"""2
Initialize arrays of walkers and their statuses3
4
N - number of agents per cluster5
centers - Center points of clusters, Tuple of tuples6
L - Box size for simulation7
"""8
function initialize(N, L, rng)9
positions = rand(rng, 1-L:L-1, (2, N)) # Spread agents out across board10
positions = positions .+ (rand(rng, 2, N) .- 0.5) # Give each agent a unique pos11
12
states = [S for i = 1:N]13
states[1] = I14
return positions, states15
end;xxxxxxxxxx19
1
"""2
Walk a single point by a random amount.3
"""4
function walk(x, y, L, rng)5
i = randn(rng)6
j = randn(rng)7
8
# move walkers by random amount9
x += i10
y += j11
12
# if the posision exceeds the boundary L, then move back13
if L ≤ abs(x) || L ≤ abs(y)14
x -= i15
y -= j16
end17
18
return x, y19
end ;xxxxxxxxxx4
1
"""2
Walk method that accepts a 1x2 array or a 2-tuple as input.3
"""4
walk(array::Union{Array, Tuple}, L, rng) = walk(array[1], array[2], L, rng);xxxxxxxxxx14
1
"""2
Get the trajectory of a random walker over N time steps.3
"""4
function trajectory(N, L=200, rng=MersenneTwister(convert(Int, 1e9*time())))5
# Pre-allocate array for efficiency6
path = Array{Tuple{Float64, Float64}, 1}(undef, N)7
# Start at the origin8
path[1] = (0.0, 0.0)9
# Do the walk10
for i = 2:N11
path[i] = walk(path[i-1], L, rng)12
end13
return path14
end;0.0
0.0
-0.568053
1.06449
-2.70724
1.03779
-3.17774
1.50317
-4.34155
1.72163
-3.72124
1.81677
-4.4771
1.59397
-4.32362
1.37935
-3.76659
2.92995
-5.62019
3.41001
-4.64491
1.92431
-4.81016
2.00098
-4.60605
0.195004
-3.43787
-0.0139538
-1.83466
1.6318
-0.661972
1.83245
-0.930706
1.56163
-2.15212
1.86684
-1.23116
2.86347
-1.11357
1.71749
-34.7818
-44.0244
-34.5511
-45.638
-35.0571
-45.5881
-35.0191
-47.1711
-34.2387
-46.5486
-35.3776
-46.6002
-34.4667
-48.6668
-34.5051
-48.4443
-35.557
-50.0915
-34.3304
-47.5218
0.0
0.0
-1.07022
0.912116
1.06969
2.54703
0.39269
2.52058
-0.00220888
3.27191
1.12974
1.42972
1.68576
2.53508
2.33011
2.87971
2.90159
2.7298
3.64599
1.60826
2.47991
2.59611
1.23927
0.178774
1.1836
0.669461
1.73043
-0.0062647
2.23047
-0.0931176
2.09465
-1.34839
2.41361
-1.58644
2.11697
-0.637231
2.24632
-0.35955
2.01647
-0.684992
-105.368
-25.6671
-104.231
-26.048
-105.569
-27.4512
-106.074
-28.2481
-105.939
-26.1177
-105.306
-26.4104
-105.003
-26.0416
-105.561
-25.8663
-105.129
-25.8418
-105.874
-25.442
0.0
0.0
0.0362361
0.649692
0.22077
-0.815234
-0.285447
-0.443169
0.210603
-1.75638
-1.08125
-2.5775
0.262841
-1.49792
-0.0602894
-2.48519
-0.217311
-1.99036
-1.56957
-2.9113
-1.20657
-2.73786
-1.44319
-1.30756
-2.91143
-2.12824
-3.83488
-1.18346
-3.57881
-2.72577
-3.31468
-3.01831
-3.82665
-3.46878
-3.23942
-3.20256
-1.00179
-3.25405
0.712356
-3.74669
98.1897
99.817
98.453
100.112
97.0904
101.267
95.6143
98.7808
96.1519
99.6258
93.6693
99.1283
94.423
98.5957
93.5126
97.4725
93.5564
96.2852
92.9709
96.314
0.0
0.0
-1.12205
0.0912732
1.698
2.0765
0.911603
2.78858
-0.521353
4.19143
-0.0963567
5.76732
-0.0146754
8.18747
-0.267524
8.9834
-1.51072
8.85686
-1.79238
9.36721
-2.14834
7.9986
-1.57924
8.41703
-1.96449
7.41184
-2.29754
7.58376
-0.656424
6.39706
-0.928671
4.62918
-1.44802
5.37714
-1.19276
7.00668
-3.05809
6.27786
-2.18688
6.34378
74.6508
-197.264
74.3404
-197.024
74.8096
-196.504
76.3885
-197.849
76.748
-197.802
75.5921
-197.981
74.519
-197.571
73.0904
-197.194
72.958
-195.755
72.4205
-196.431
0.0
0.0
1.58481
-1.19239
1.09027
-0.902452
0.729667
-0.259724
-0.0472274
1.17939
0.781782
1.91744
-0.616017
4.44143
0.95741
4.98606
0.600343
4.61267
-0.388462
6.77372
0.113766
6.21515
-0.433426
5.38151
-1.96368
5.44272
-2.63772
4.9752
-2.52101
6.68561
-5.11613
8.90105
-4.46547
8.60541
-4.45178
8.98679
-5.71487
8.78463
-6.19819
9.21916
22.6555
160.79
23.0811
161.525
22.4738
160.287
24.875
159.248
25.6051
160.724
24.9761
160.648
23.5865
161.145
21.971
160.986
22.488
161.885
22.4092
160.623
0.0
0.0
0.153276
-0.660846
-0.112499
-1.62266
-0.282843
-1.19401
0.873319
-0.967179
-0.606692
-1.69308
-0.875057
0.0480284
-0.798583
-0.230233
-1.03216
-0.865882
-1.15041
-0.477491
-1.7122
-0.511556
0.566429
-0.412889
-0.135687
-1.22445
-1.17716
0.367942
-0.248955
1.17221
-0.358074
1.15429
-0.432614
0.752704
0.145676
0.515431
0.923263
-1.80629
1.26642
-3.25731
92.15
-72.3255
93.8048
-73.6655
93.5581
-75.9972
93.7374
-78.2687
93.2459
-78.9634
93.3564
-77.0435
91.9846
-77.2658
95.2237
-76.908
97.5896
-75.317
97.8566
-76.4158
0.0
0.0
0.327544
-0.930099
0.401632
0.435788
0.925549
0.575476
0.256502
-0.355065
-0.322153
-0.79379
0.705694
-0.868012
0.890335
-1.29825
0.422843
-2.47183
0.163403
-1.65528
-1.40554
-1.44282
-0.579072
-0.38709
1.40706
1.3466
2.09734
0.881204
1.74244
-1.42896
0.838285
-0.320601
1.4749
-1.32724
1.79429
-0.0112178
3.81011
0.0568438
2.90132
0.0153873
12.6133
-50.8934
13.3764
-50.6184
13.8414
-50.6385
14.4069
-51.4919
14.1699
-51.3052
13.2196
-50.5674
13.9206
-50.3139
15.1482
-49.8618
13.6588
-49.3243
14.1506
-50.0092
0.0
0.0
-0.664554
0.597482
1.10206
3.17821
1.63944
2.38021
0.878587
2.49472
2.27361
3.19458
2.60113
2.36175
3.07954
3.24994
3.11898
3.50453
4.05079
4.09663
5.19058
4.37908
5.98779
2.80146
5.47834
3.78884
4.71282
2.93023
5.74461
4.92447
5.4782
2.34982
5.69938
2.29695
7.67557
2.64349
9.01142
3.95696
9.41019
5.21716
-29.5425
107.706
-30.3915
107.181
-30.7365
106.089
-30.0148
107.445
-31.3888
106.439
-30.9523
106.205
-28.8105
106.792
-27.1216
106.543
-27.7227
106.728
-28.807
104.743
0.0
0.0
-0.614896
2.67956
0.0199324
2.97126
0.0344443
4.66061
-0.253578
3.44471
-0.692148
4.08203
-0.746913
4.97167
0.547313
4.21671
-0.0405518
4.88547
-1.06908
3.68432
-0.798227
3.47161
0.112525
3.80926
0.766659
3.92909
-0.252964
4.33525
0.910638
4.51406
0.577749
3.39757
0.545448
4.02388
-1.45918
2.85628
-3.12704
0.839895
-2.43633
0.43517
111.713
-44.568
111.209
-46.7496
110.935
-47.6315
111.43
-48.8843
112.046
-47.257
111.744
-47.6466
113.092
-47.1131
113.008
-45.2413
112.591
-45.0062
112.419
-45.2611
0.0
0.0
-0.66536
-0.887734
-0.242284
-0.464673
-2.12131
-2.56266
-2.49762
-1.38261
-2.43041
-0.105131
-0.335501
0.161143
-0.632506
0.818967
-0.813012
0.154655
-0.37407
-1.16996
-0.0994359
-2.70629
-1.03473
-1.36218
-0.148491
-2.6741
-0.158258
-0.454979
1.09792
0.241442
1.20942
1.35246
1.24796
0.582064
2.50575
-2.3604
2.01385
-0.540331
2.89252
-1.12238
54.5466
130.991
54.5425
131.364
54.0648
131.635
54.4371
130.511
57.4385
131.086
57.7913
132.684
58.7007
133.029
59.4782
130.814
57.1397
130.837
55.5555
132.721
xxxxxxxxxx2
1
# Make an array of 10 walkers2
walkers = [trajectory(20_000) for i in 1:10]The Simulation
Define close contacts using NearestNeighbors.jl
xxxxxxxxxx4
1
"""2
Convenience function combining `BallTee` creation and `inrange` functions.3
"""4
encounters(contacts, individual, radius) = inrange(BallTree(contacts), individual, radius);step! sweep! sim!
One step! moves one agent in one time step.
One sweep! moves all agents on average in one time step
One simulation is a series of time sweeps
xxxxxxxxxx34
1
function step!(positions, states, N, L, pᵢ, pᵣ, radius, rng::MersenneTwister)2
# Pick an agent at random3
i = rand(rng, 1:N)4
# Give infected agents a chance to recover before infecting others5
if states[i] == I6
if rand(rng) ≤ pᵣ7
states[i] = R8
end9
end10
11
# Move the agent to a new location12
xy = walk(positions[:, i], L, rng)13
positions[1, i] = xy[1]14
positions[2, i] = xy[2]15
16
# If the agent is still infected17
if states[i] == I18
# Determine close contacts19
contacts = encounters(positions, positions[:, i], radius)20
# Remove self21
contacts = contacts[(contacts .≠ i)]22
# For each contact, 23
for contact in contacts24
# Give a chance to infect if susceptible25
if states[contact] == S26
if rand(rng) ≤ pᵢ27
states[contact] = I28
end29
end30
end31
end32
33
return positions, states34
end;xxxxxxxxxx6
1
function sweep!(positions, states, N, L, pᵢ, pᵣ, radius, rng::MersenneTwister)2
for n = 1:N3
step!(positions, states, N, L, pᵢ, pᵣ, radius, rng)4
end5
return positions, states6
end;xxxxxxxxxx37
1
function simulation(Sw, N, L, pᵢ, pᵣ, radius, rng=MersenneTwister(convert(Int, 1e9*time())))2
system = Dict(3
"S" => Sw, # Number of sweeps ("days")4
"N" => N, # Number of agents5
"L" => L, # Length of bounding box for agents6
"pᵢ" => pᵢ, # Infection probability7
"pᵣ" => pᵣ, # Recovery probability8
"positions" => [], # Positions per sweep9
"states" => [], # States per sweep10
"susceptible" => Int[], # Number of susceptible agents at step n11
"infected" => Int[], # Number of infected agents at step n12
"recovered" => Int[]) # Number of recovered agents at step n13
14
15
positions, states = initialize(N, L, rng)16
# Large sim threshold for saving movement history17
thresh = 1e518
if Sw*N < thresh19
push!(system["positions"], deepcopy(positions))20
end21
push!(system["states"], deepcopy(states))22
push!(system["susceptible"], sum(states .== S))23
push!(system["infected"], sum(states .== I))24
push!(system["recovered"], sum(states .== R))25
26
for n=2:Sw27
sweep!(positions, states, N, L, pᵢ, pᵣ, radius, rng)28
if Sw*N < thresh29
push!(system["positions"], deepcopy(positions))30
end31
push!(system["states"], deepcopy(states))32
push!(system["susceptible"], sum(states .== S))33
push!(system["infected"], sum(states .== I))34
push!(system["recovered"], sum(states .== R))35
end36
return system37
end;Visualize it
xxxxxxxxxx1
1
med_sim = simulation(100, 800, 20, 0.3, 0.05, 1);xxxxxxxxxx9
1
# "Big" Simulation parameters2
begin3
days = 5004
N = 15005
boxsize = 30 6
prob_i = 0.17
prob_r = 0.028
cluster_rad = 2.59
end;"susceptible"
1499
1498
1497
1495
1493
1491
1483
1471
1448
1420
1401
1376
1356
1331
1316
1297
1272
1240
1206
1166
0
0
0
0
0
0
0
0
0
0
"recovered"
0
0
0
0
0
0
0
1
1
3
6
9
13
15
18
22
25
31
32
38
1500
1500
1500
1500
1500
1500
1500
1500
1500
1500
"S"
500
"infected"
1
2
3
5
7
9
17
28
51
77
93
115
131
154
166
181
203
229
262
296
0
0
0
0
0
0
0
0
0
0
"N"
1500
"pᵣ"
0.02
"L"
30
"states"
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
R::InfectionStatus = 2
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
I::InfectionStatus = 1
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
S::InfectionStatus = 0
I::InfectionStatus = 1
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
R::InfectionStatus = 2
"positions"
"pᵢ"
0.1
xxxxxxxxxx1
1
sim = simulation(days, N, boxsize, prob_i, prob_r, cluster_rad)The SIR Model
Has three states:
Susceptible
Infectious
Recovered
The time evolution of these states can be described by a set of ordinary differential equations.
The approximate solution can be taken in small time steps.
xxxxxxxxxx25
1
"""2
Get the approximate solution to the SIR ODE equations.3
"""4
function euler_SIR(S₀, I₀, R₀, β, γ, h, T)5
# get the total number of timesteps to use6
n = ceil(Int, T/h)7
# Pre-allocate the arrays8
t = collect(1:h:n)9
10
S = Array{Float64, 1}(undef, n)11
I = Vector{Float64}(undef, n)12
R = zeros(Float64, n)13
14
S[1] = S₀15
I[1] = I₀16
R[1] = R₀17
18
for i = 1:(n-1)19
t[i+1] = t[i] + h20
S[i+1] = S[i] - h*(β*S[i]*I[i])21
I[i+1] = I[i] + h*(β*S[i]*I[i] - γ*I[i])22
R[i+1] = R[i] + h*γ*I[i]23
end24
return [t, S, I, R]25
end;1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0
13.0
14.0
15.0
16.0
17.0
18.0
19.0
20.0
491.0
492.0
493.0
494.0
495.0
496.0
497.0
498.0
499.0
500.0
0.9999
0.99988
0.999856
0.999828
0.999794
0.999754
0.999706
0.99965
0.999582
0.999502
0.999406
0.999292
0.999157
0.998996
0.998804
0.998576
0.998305
0.997983
0.997599
0.997143
9.87317e-10
9.84874e-10
9.82462e-10
9.8008e-10
9.77728e-10
9.75405e-10
9.7311e-10
9.70844e-10
9.68605e-10
9.66394e-10
0.0001
0.000118998
0.000141605
0.000168506
0.000200516
0.000238606
0.000283929
0.000337859
0.000402028
0.00047838
0.000569225
0.00067731
0.000805903
0.000958888
0.00114088
0.00135738
0.0016149
0.00192118
0.00228543
0.00271856
0.0123686
0.0122449
0.0121225
0.0120012
0.0118812
0.0117624
0.0116448
0.0115283
0.011413
0.0112989
0.0
1.0e-6
2.18998e-6
3.60603e-6
5.29108e-6
7.29624e-6
9.6823e-6
1.25216e-5
1.59002e-5
1.99205e-5
2.47043e-5
3.03965e-5
3.71696e-5
4.52286e-5
5.48175e-5
6.62264e-5
7.98001e-5
9.59491e-5
0.000115161
0.000138015
0.987631
0.987755
0.987878
0.987999
0.988119
0.988238
0.988355
0.988472
0.988587
0.988701
xxxxxxxxxx1
1
test = euler_SIR(0.9999, 0.0001, 0.0, 0.2, 0.01, 1, 500)Optimizing with Gradient Descent
Use derivatives of a function to find its mininmum.
The derivative of a function can be approximated as:
f (generic function with 1 method)xxxxxxxxxx1
1
f(x) = x^3 - x^2 + x - 1deriv (generic function with 2 methods)xxxxxxxxxx1
1
deriv(f, a, h=1.0e-5) = (f(a+h) - f(a))/h # Forward differencetangent_line (generic function with 1 method)xxxxxxxxxx8
1
function tangent_line(f, x, a; h=1.0e-5)2
# y = m*(x-a) + b(a) <- equation of line through a3
df = x -> deriv(f, a, h) # for all values of x, input deriv at `a`4
b = f(a) # y-intercept at a5
m = df(a) # slope at a6
7
return m*(x-a) + b # Equation of line at a8
enda =
Extending to 2-D with the gradient, ∇
For a function g(x, y), the gradient is a vector of partial derivatives with respect to x and y.
xxxxxxxxxx1
1
g(x, y) = x^2 + y^3;xxxxxxxxxx1
1
dx(f, a, b; h= 1e-5) = deriv(x -> f(x, b), a, h);xxxxxxxxxx1
1
dy(f, a, b; h= 1e-5) = deriv(y -> f(a, y), b, h);xxxxxxxxxx1
1
gradient(f, a, b; h=1e-5) = [dx(f, a, b, h=h), dy(f, a, b, h=h)];4.00001
75.0001
xxxxxxxxxx1
1
gradient(g, 2, 5)Gradient descent optimization in 1-D
xxxxxxxxxx33
1
function descent_1d(f, x0, η=0.01; max_iter=100_000)2
# Calculate the derivative of the input function3
df = y -> deriv(f, y)4
5
# x-values of the descent6
x = Float64[x0]7
8
count = 1 # Set a counter to help break the loop9
while count < max_iter10
# Chaning the sign of η simplifies the code later11
ξ = (0 < df(x[end])) ? -η : η12
13
# Trying to get around inflection point problem14
# by computing ahead in batches15
xᵢ = [x[end] + n*ξ for n = 1:5] # next 10 points16
fxᵢ = f.(xᵢ) # Calculate the function value for the points17
dfdxᵢ = abs.(df.(xᵢ)) # Calculate the derivatives for the points18
19
# Find the index of the minimum value of the derivative20
min_ind = findfirst(dfdxᵢ .== minimum(dfdxᵢ)) # Min dfdx index21
22
# Add x values to the array up to the minimum df/dx value23
x = vcat(x, xᵢ[1:min_ind])24
25
# If the smallest value of f(x) is not at the end,26
# then the minimum has already been reached.27
fxᵢ[end] ≠ minimum(fxᵢ) && break28
29
count += 130
end31
32
return x33
end;Let's try it out!
xxxxxxxxxx1
1
foo(x) = x^4 + 3x^3 + 5;1.0
0.99
0.98
0.97
0.96
0.95
0.94
0.93
0.92
0.91
0.9
0.89
0.88
0.87
0.86
0.85
0.84
0.83
0.82
0.81
-2.16
-2.17
-2.18
-2.19
-2.2
-2.21
-2.22
-2.23
-2.24
-2.25
xxxxxxxxxx1
1
x_descent = descent_1d(foo, x0)9.0
8.87149
8.74594
8.62331
8.50355
8.38663
8.2725
8.16112
8.05246
7.94646
7.8431
7.74233
7.64411
7.54841
7.45518
7.36438
7.27598
7.18994
7.10623
7.02479
-3.46526
-3.4812
-3.49539
-3.5078
-3.5184
-3.52715
-3.53402
-3.53897
-3.54196
-3.54297
xxxxxxxxxx1
1
y_descent = foo.(x_descent)x₀ =
index = 1
Carl Friedrich Gauss walks into a bar...
and tries to fit his data.
xxxxxxxxxx1
1
gaussian(x, μ=0, σ=1) = 50/(σ*√(2π))*ℯ^(-(x - μ)^2/(2σ^2));xxxxxxxxxx1
1
noisy_gaussian(x, μ=0, σ=1) = gaussian(x, μ, σ) + randn();xxxxxxxxxx1
1
histogram(randn(100_000))Gradient descent optimization in 2-D
xxxxxxxxxx38
1
function gradient_descent_2d(f, x0, y0; n=100_000)2
# Define the step size3
η = 1.0e-44
# Pre-allocate an array of ones (will be overwritten)5
descent = ones(Float64, (n, 2))6
# Add the starting point to the first array entry7
descent[1, :] = [x0, y0]8
# This step avoids cyclic references (a Pluto thing)9
max_i = copy(n)10
# Loop through until the minimum is found11
# OR max iterations is reached12
for i=1:n-113
# Calculate the gradient of the point14
grad = gradient(f, descent[i, :]...)15
# For each direction (x, y)16
for j=1:217
if grad[j] == 0 # If the grad is zero, 18
descent[i+1, j] = descent[i, j] # Just add the point19
elseif grad[j] > 0 # If the grad is +, 20
descent[i+1, j] = descent[i, j] - η # Then step backward21
elseif grad[j] < 0 # If the grad is -, 22
descent[i+1, j] = descent[i, j] + η # Then step forward23
end24
end25
26
if i > 227
# If the point tries to oscillate back-and-forth,28
# then it's a minimum, so quit.29
pos1 = descent[i, :] # Current position30
pos2 = descent[i-2, :] # Position two steps ago31
if pos1 == pos2 # If they're equal,32
max_i = i-2 # then the minimum has been reached33
break # so quit.34
end35
end36
end37
return descent[1:max_i, :]38
end;Use gradient descent to "learn" μ and σ
xxxxxxxxxx1
1
mu = 5; sigma = 0.6;-10.0:0.005:10.0xxxxxxxxxx1
1
xg = -10:0.005:106.37827e-135
7.85536e-135
9.67384e-135
1.19125e-134
1.46681e-134
1.80599e-134
2.22346e-134
2.73722e-134
3.36948e-134
4.14748e-134
5.10476e-134
6.28256e-134
7.73157e-134
9.51413e-134
1.17068e-133
1.44039e-133
1.77211e-133
2.18008e-133
2.68177e-133
3.29869e-133
5.15551e-14
4.81248e-14
4.49196e-14
4.19249e-14
3.91272e-14
3.65137e-14
3.40723e-14
3.1792e-14
2.96622e-14
2.76732e-14
xxxxxxxxxx1
1
yg = gaussian.(xg, mu, sigma)-0.694508
-1.44066
0.0499991
0.273882
14.5429
31.7296
12.0714
-1.13436
0.146313
-0.649156
0.176367
xxxxxxxxxx4
1
begin2
xng = sort(collect(0:1:10) .+ (rand(11)*2 .- 1))3
yng = noisy_gaussian.(xng, mu, sigma)4
endThe loss function measures how much the model deviates from the data.
xxxxxxxxxx1
1
loss(μ=0, σ=1) = log(sum((gaussian.(xng, μ, σ) .- yng).^2)); # Log of sum squared errorxxxxxxxxxx1
1
xμ = -2:0.1:10;xxxxxxxxxx1
1
yσ = 0.2:0.1:6.0;z (generic function with 1 method)xxxxxxxxxx1
1
z(b, a) = loss(a, b)Passing loss into gradient_descent selects parameters that drive it towards zero.
xxxxxxxxxx1
1
params = gradient_descent_2d(loss, 0, 1);step =
How does our SIR model measure up?
Need to redo the loss function
xxxxxxxxxx1
1
sse(x_data, x_theory) = sum((x_theory .- x_data).^2); # Sum of squared errorxxxxxxxxxx10
1
function loss_sir(β, γ) 2
model_values = euler_SIR(0.9999, 0.0001, 0.0, β, γ, 1, length(sim["susceptible"]))3
4
# One loss value for each state5
Lₛ = sse(sim["susceptible"], model_values[2])6
Lᵢ = sse(sim["infected"], model_values[3])7
Lᵣ = sse(sim["recovered"], model_values[4])8
9
return Lₛ + 4*Lᵢ + Lᵣ # Put the most weight on the infectious curve10
end;3186×2 Matrix{Float64}:
0.1 0.0
0.1001 0.0001
0.1002 0.0002
0.1003 0.0003
0.1004 0.0004
0.1005 0.0005
0.1006 0.0006
⋮
0.418 0.0164
0.4181 0.0163
0.4182 0.0164
0.4183 0.0163
0.4184 0.0164
0.4185 0.0163xxxxxxxxxx1
1
ode_params = gradient_descent_2d(loss_sir, beta, gamma, n=10_000)Fit Parameters
β = 0.4185
γ = 0.0163
step = 1
β =